      ! Returns the next concentration in the array (taking boundary conditions into
      ! account)
      subroutine get_next_val(next_x, x, i, n)
          implicit none;
          INTEGER, INTENT(in) :: i, n;
          DOUBLE PRECISION, INTENT(inout) :: next_x;
          DOUBLE PRECISION, INTENT(in), DIMENSION(n) :: x;
          IF (i == 1) THEN
              next_x = x(i+1);
          ELSE IF (i == n) THEN
              next_x = 1.0; ! Constant outflux
          ELSE
              next_x = x(i+1);
          END IF
      END SUBROUTINE get_next_val

